// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtMOABFileFormat.C                           //
// ************************************************************************* //

#include <avtMOABFileFormat.h>

#include <string>

#include <vtkCellType.h>
#include <vtkFloatArray.h>
#include <vtkIntArray.h>
#include <vtkBitArray.h>

#include <vtkUnstructuredGrid.h>

#include <avtCallback.h>
#include <avtDatabaseMetaData.h>

#include <DBOptionsAttributes.h>
#include <DebugStream.h>
#include <Expression.h>

#include <InvalidVariableException.h>

// MOAB includes
#include <moab/Core.hpp>
#include <moab/EntityType.hpp>

#include <moab/mhdf_public.h>
#include <moab/VtkUtil.hpp>

// parallel stuff
#ifdef PARALLEL
#include <mpi.h>
#include <avtParallel.h>
#include <moab/ParallelComm.hpp>
#endif
using     std::string;

//using namespace moab;

#define MBVIS_CHK_ERR(err_code) \
  do { \
    if (moab::MB_SUCCESS != err_code) { \
        std::cout << "My error code is " << err_code << " and I failed this plugin" << std::endl; \
        moab::MBError(__LINE__, __func__, __FILENAME__, __MBSDIR__, moab::MB_FAILURE, "MOAB reader failure.", moab::MB_ERROR_TYPE_NEW_LOCAL); \
        throw err_code; \
    } \
  } while(false)


// ****************************************************************************
//  Method: avtMOABFileFormat constructor
//
//  Programmer: vijaysm -- generated by xml2avt
//  Creation:   Wed Jan 20 13:02:35 PST 2016
//
// ****************************************************************************

avtMOABFileFormat::avtMOABFileFormat(const char *filename, const DBOptionsAttributes *readOpts)
    : avtSTMDFileFormat(&filename, 1), readOptions(readOpts),  file_descriptor(NULL), pcomm(NULL), mbCore(NULL)
{
    // INITIALIZE DATA MEMBERS

	debug1<< " constructor called, file " << filename << "\n";
    fileName = strdup(filename);
    opt1d = false;
    opt2d = true;
    opt3d = true;
    edges = new moab::Range;
    faces = new moab::Range;
    solids = new moab::Range;
    select = new moab::Range;
}


// ****************************************************************************
//  Method: avtMOABFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: vijaysm -- generated by xml2avt
//  Creation:   Wed Jan 20 13:02:35 PST 2016
//
// ****************************************************************************

avtMOABFileFormat::~avtMOABFileFormat() {
    debug1 << " avtMOABFileFormat::~avtMOABFileFormat for file:" << fileName << "\n";
    FreeUpResources();
    if (fileName != NULL)
    {
        free(fileName);
        fileName = NULL;
    }
}
void
avtMOABFileFormat::FreeUpResources(void) {
    debug1 << " avtMOABFileFormat::FreeUpResources: freeing file descriptor for file " << fileName << "\n";
    if (file_descriptor) {
        free(file_descriptor);
        file_descriptor = NULL;
    }
#ifdef PARALLEL
	if (pcomm)
	{
		delete pcomm;
		pcomm = NULL;
	}
#endif
    if (mbCore) {
        delete edges;
        delete faces;
        delete solids;
        delete select;
        delete mbCore;
        mbCore = NULL;
    }
}

void
avtMOABFileFormat::gatherMhdfInformation()
{
    //MHDF_FileDesc* file_descriptor;
    unsigned long size = 0;
    if (0==rank)
    {
        MHDF_FileHandle file;
        MHDF_Status status;
        unsigned long max_id;

        file = mhdf_openFile( fileName, 0, &max_id, -1, &status );
        if (mhdf_isError( &status )) {
            debug1 << "fail to open file " << fileName << "\n";
        }
        else
        {
            debug1 << "opened file " << fileName << "\n";
        }
        file_descriptor = mhdf_getFileSummary( file, H5T_NATIVE_ULONG, &status, 1);
        debug1 << " avtMOABFileFormat::gatherMhdfInformation file_descriptor address: " << (void*) file_descriptor << "\n";
        if (mhdf_isError( &status )) {
            debug1 << "fail to get summary\n";
        }
        else
        {
            debug1 << "got summary\n";
        }
        size = file_descriptor->total_size;
        debug1 << "file descriptor size =" << size << "\n";
        file_descriptor->offset = (unsigned char*)file_descriptor;
        mhdf_closeFile( file, &status );
    }
#ifdef PARALLEL
    // need to communicate file summary to all processes
    int mpi_err = MPI_Bcast(&size, 1, MPI_UNSIGNED_LONG, 0, VISIT_MPI_COMM);
    if (mpi_err || !size)
    {
        debug1 << "fail to communicate size: " << size << "\n";
    }
    if (0 != rank)
        file_descriptor = reinterpret_cast<MHDF_FileDesc*>(malloc(size));
    mpi_err = MPI_Bcast(file_descriptor, size, MPI_BYTE, 0, VISIT_MPI_COMM);
    if (mpi_err )
    {
        debug1 << "fail to communicate file summary: \n";
    }

    if (0 != rank)
        mhdf_fixFileDesc(file_descriptor, reinterpret_cast<MHDF_FileDesc*>(file_descriptor->offset));
#endif
    // extract number of nodes, elements, dense tags for the time being
    long num_nodes = file_descriptor->nodes.count;
    int number_node_tags = file_descriptor->nodes.num_dense_tags;

    debug2 << "Nodes: " << num_nodes << " dense tags: " <<  number_node_tags << "\n";
    for (int i=0; i< number_node_tags; i++)
    {
        MHDF_TagDesc & tagStr = file_descriptor->tags[file_descriptor->nodes.dense_tag_indices[i]];
        const char * tag_name = tagStr.name;

        int sizeTag = tagStr.size;
        debug2 << "   tag "<< i << " "  << tag_name << " size:" << sizeTag << "\n";
        struct tagBasic  tag1;
        tag1.nameTag=string(tag_name);
        tag1.size=sizeTag;
        if (tagStr.default_value && ( sizeTag == 1 ) &&
                ( (tagStr.type == mhdf_INTEGER) || (tagStr.type == mhdf_FLOAT )) )
        {
            tag1.defValue = tagStr.default_value;
            tag1.type = tagStr.type;
        }

        debug2 << "  node     tag "<< i << " "  << tag_name <<" size:" << sizeTag << "\n";
        nodeTags.push_back(tag1);
    }

    debug2 << "Cells descriptors " << file_descriptor->num_elem_desc <<  "\n";

    for (int i=0; i< file_descriptor->num_elem_desc; i++)
    {
        const char * etype = file_descriptor->elems[i].type;

        MHDF_EntDesc & edesc = file_descriptor->elems[i].desc;
        int number_elem_tags = edesc.num_dense_tags;
        debug2 << "   elem type  "<< i << " "  << etype << " num dense tags: "
            << number_elem_tags <<"\n";
        for (int j=0; j<number_elem_tags; j++)
        {
            MHDF_TagDesc & tagStr = file_descriptor->tags[edesc.dense_tag_indices[j]];
            const char * tag_name = tagStr.name;
            int sizeTag = tagStr.size;
            struct tagBasic  tag1;

            tag1.nameTag=string(tag_name);
            tag1.size=sizeTag;
            if (tagStr.default_value && ( sizeTag == 1 ) &&
                            ( (tagStr.type == mhdf_INTEGER) || (tagStr.type == mhdf_FLOAT )) )
            {
                tag1.defValue = tagStr.default_value;
                tag1.type = tagStr.type;
            }

            debug2 << " elem   tag "<< j << " "  << tag_name <<" size:" << sizeTag << "\n";
            elemTags.insert(tag1);
        }
    }

    num_parts = file_descriptor->numEntSets[0];
    num_mats  = file_descriptor->numEntSets[1];
    num_neumann = file_descriptor->numEntSets[2];
    num_diri = file_descriptor->numEntSets[3];
    num_geom = file_descriptor->numEntSets[4];

    // gather material values ids in the set
    for (int i=0; i<num_mats; i++)
      materials.insert(file_descriptor->defTagsVals[1][i]);

    if (nProcs > 1 && num_parts < nProcs)
    {
        // EXCEPTION1(InvalidVariableException, "too few parts in PARALLEL_PARTITION ");
        debug1 << " there are less parts in partition than processors " << "nProcs=" << nProcs << " num_parts:" << num_parts << " \n";
    }

    return;
}

// ****************************************************************************
//  Method: avtMOABFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: vijaysm -- generated by xml2avt
//  Creation:   Wed Jan 20 13:02:35 PST 2016
//
// ****************************************************************************

void
avtMOABFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
{
    rank=0;
    nProcs=1;
#ifdef PARALLEL
    MPI_Comm_rank(VISIT_MPI_COMM, &rank);
    MPI_Comm_size(VISIT_MPI_COMM, &nProcs);
#endif
    debug1 << "avtMOABFileFormat::PopulateDatabaseMetaData rank:" << rank << "\n";

    try
    {

        // we find a lot of info about the mesh, without actually loading it
        gatherMhdfInformation();

        moab::ErrorCode merr;
        std::string domainSetClassName;

        domainSetClassName = "Root";

        //
        // Define the mesh
        //
        string meshname = "mesh";
        avtMeshMetaData *mesh = new avtMeshMetaData;
        mesh->name = meshname;
        mesh->meshType = AVT_UNSTRUCTURED_MESH;
        // we are calling later md->SetFormatCanDoDomainDecomposition(true)
        //  because of that,
        mesh->numBlocks = 1; // nProcs;
        mesh->blockOrigin = 1;
        mesh->spatialDimension = file_descriptor->nodes.vals_per_ent; // usually 3, but some files might have 2
        mesh->topologicalDimension = 3; // we need to look at element types
        mesh->blockTitle = "Blocks";
        mesh->blockPieceName = domainSetClassName; // "Block%012d";
        mesh->hasSpatialExtents = false;

        //
        // Hack to provide feedback we're running correct plugin
        //
        md->SetDatabaseComment("This is the New MOAB Plugin");

        if (md != NULL)
            md->Add(mesh);

        for (std::set<struct tagBasic>::iterator setIter = elemTags.begin(); setIter!=elemTags.end(); setIter++)
        {
            string nameDisplay("ELEM_");
            struct tagBasic tag = *setIter;
            nameDisplay += tag.nameTag;
            if (tag.size == 1)
            {
                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = nameDisplay.c_str();
                smd->meshName = meshname;
                smd->centering = AVT_ZONECENT;
                double defValDouble = 0;
                if (tag.defValue && tag.type == 1) // integer
                {
                    defValDouble = *((int*)(tag.defValue));
                }
                if (tag.defValue && tag.type == 2) // double
                {
                    defValDouble = *((double*)(tag.defValue));
                }
                double values[2] = {defValDouble, defValDouble};
                smd->SetMissingData(values);
                smd->SetMissingDataType(avtScalarMetaData::MissingData_Value);
                md->Add(smd);
            }
            else if (tag.size == 2 || tag.size == 3)
              AddVectorVarToMetaData(md, nameDisplay.c_str(), meshname, AVT_ZONECENT, tag.size);
            else
              AddArrayVarToMetaData(md, nameDisplay.c_str(), tag.size, meshname, AVT_ZONECENT);
        }
        for (int i=0; i<nodeTags.size(); i++)
        {
            string nameDisplay("NODE_");
            struct tagBasic tag = nodeTags[i];
            nameDisplay +=tag.nameTag;
            if (tag.size == 1)
            {
                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = nameDisplay.c_str();
                smd->meshName = meshname;
                smd->centering = AVT_NODECENT;
                double defValDouble = 0;
                if (tag.defValue && tag.type == 1) // integer
                {
                    defValDouble = *((int*)(tag.defValue));
                }
                if (tag.defValue && tag.type == 2) // double
                {
                    defValDouble = *((double*)(tag.defValue));
                }
                double values[2] = {defValDouble, defValDouble};
                smd->SetMissingData(values);
                smd->SetMissingDataType(avtScalarMetaData::MissingData_Value);
                md->Add(smd);
            }
              // AddScalarVarToMetaData(md, nameDisplay.c_str(), meshname, AVT_NODECENT);
            else if (tag.size == 2 || tag.size == 3)
              AddVectorVarToMetaData(md, nameDisplay.c_str(), meshname, AVT_NODECENT, tag.size);
            else
              AddArrayVarToMetaData(md, nameDisplay.c_str(), tag.size, meshname, AVT_NODECENT);
        }

        //  So, here, we handle the parallel partition
        // decomposition as an enumerated scalar variable.
        //  assume that parallel partitions numbers start at 0 and end at num_parts -1
        if (num_parts>0)
        {
            avtScalarMetaData *ppsmd = new avtScalarMetaData("ParallelPartition", meshname, AVT_ZONECENT);
            ppsmd->SetEnumerationType(avtScalarMetaData::ByValue);
            for (int i = 0 ; i < num_parts ; i++)
            {
                char name[256];

                snprintf(name, sizeof(name), "PAR_PART_%d", i);
                ppsmd->AddEnumNameValue(name, i+1);
            }
            ppsmd->AddEnumNameValue("NOT_SPECIFIED", num_parts+1);
            md->Add(ppsmd);
        }
        if (materials.size()>0)
        {
            avtScalarMetaData *msmd = new avtScalarMetaData("Materials", meshname, AVT_ZONECENT);
            msmd->SetEnumerationType(avtScalarMetaData::ByValue);
            int i=1;
            msmd->AddEnumNameValue("NOT_SPECIFIED", 0);
            for (std::set<int>::iterator mIt = materials.begin();
                mIt!=materials.end(); mIt++, i++)
            {
                int value= *mIt;
                char name[256];
                snprintf(name, sizeof(name), "Material_%d", value);
                msmd->AddEnumNameValue(name, value);
            }
            md->Add(msmd);
        }
        if (num_neumann>0)
        {
            avtScalarMetaData *nsmd = new avtScalarMetaData("NeumannSets", meshname, AVT_ZONECENT);
            nsmd->SetEnumerationType(avtScalarMetaData::ByBitMask);
            for(int j=0; j<num_neumann; j++)
                neumannsets.insert(file_descriptor->defTagsVals[2][j] );
            int i=1;
            for (std::set<int>::iterator mIt = neumannsets.begin();
                mIt!=neumannsets.end(); mIt++, i++)
            {
                int value= *mIt;
                char name[256];
                snprintf(name, sizeof(name), "NeumannSet_%d", value);
                nsmd->AddEnumNameValue(name, i);
            }

            // add one more universal, for all entities
            nsmd->AddEnumNameValue("NeumannUniversal", num_neumann+1);

            md->Add(nsmd);
        }
        if (num_diri>0)
        {
            avtScalarMetaData *dsmd = new avtScalarMetaData("DirichletSets", meshname, AVT_NODECENT);
            dsmd->SetEnumerationType(avtScalarMetaData::ByBitMask);
            for(int j=0; j<num_diri; j++)
                dirichsets.insert(file_descriptor->defTagsVals[3][j] );
            int i=1;
            for (std::set<int>::iterator mIt = dirichsets.begin();
                            mIt!=dirichsets.end(); mIt++, i++)
            {
                int value= *mIt;
                char name[256];
                snprintf(name, sizeof(name), "Dirichlet_%d", value);
                dsmd->AddEnumNameValue(name, i);
            }

            // add one more universal
            dsmd->AddEnumNameValue("DirichletUniversal", num_diri+1);

            md->Add(dsmd);
        }

        if (num_geom>0)
        {
            avtScalarMetaData *dsmd = new avtScalarMetaData("GeometrySets", meshname, AVT_ZONECENT);
            dsmd->SetEnumerationType(avtScalarMetaData::ByBitMask);
            for(int i=0; i<num_geom; i++)

            {
                char name[256];
                snprintf(name, sizeof(name), "GeomSet_%d_dim_%d", file_descriptor->defTagsEntSets[4][i],
                    file_descriptor->defTagsVals[4][i]);
                dsmd->AddEnumNameValue(name, i);
            }

            // add one more universal
            dsmd->AddEnumNameValue("GeometryUniversal", num_geom+1);

            md->Add(dsmd);
        }



        // Because its so convenient, add the zonetype expression as a
        // variable on the mesh
        //
        // Expression expr;
        // expr.SetName("zonetype_rank");
        // expr.SetDefinition("zonetype_rank(mesh)");
        // expr.SetType(Expression::ScalarMeshVar);
        // md->AddExpression(&expr);

        md->SetFormatCanDoDomainDecomposition(true);

    }
    catch (moab::ErrorCode errCode)
    {
        std::string errInfo;
        mbCore->get_last_error(errInfo);
        if (!avtCallback::IssueWarning(errInfo.c_str()))
            std::cerr << errInfo << std::endl;
    }

    //:q

    // CODE TO ADD A MESH
    //
    // string meshname = ...
    //
    // AVT_RECTILINEAR_MESH, AVT_CURVILINEAR_MESH, AVT_UNSTRUCTURED_MESH,
    // AVT_POINT_MESH, AVT_SURFACE_MESH, AVT_UNKNOWN_MESH
    // avtMeshType mt = AVT_RECTILINEAR_MESH;
    //
    // int nblocks = YOU_MUST_DECIDE;
    // int block_origin = 0;
    // int spatial_dimension = 2;
    // int topological_dimension = 2;
    // double *extents = NULL;
    //
    // Here's the call that tells the meta-data object that we have a mesh:
    //
    // AddMeshToMetaData(md, meshname, mt, extents, nblocks, block_origin,
    //                   spatial_dimension, topological_dimension);
    //

    //
    // CODE TO ADD A SCALAR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddScalarVarToMetaData(md, varname, mesh_for_this_var, cent);
    //

    //
    // CODE TO ADD A VECTOR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    // int vector_dim = 2;
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddVectorVarToMetaData(md, varname, mesh_for_this_var, cent,vector_dim);
    //

    //
    // CODE TO ADD A TENSOR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    // int tensor_dim = 9;
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddTensorVarToMetaData(md, varname, mesh_for_this_var, cent,tensor_dim);
    //

    //
    // CODE TO ADD A MATERIAL
    //
    // string mesh_for_mat = meshname; // ??? -- could be multiple meshes
    // string matname = ...
    // int nmats = ...;
    // vector<string> mnames;
    // for (int i = 0 ; i < nmats ; i++)
    // {
    //     char str[32];
    //     sprintf(str, "mat%d", i);
    //     -- or -- 
    //     strcpy(str, "Aluminum");
    //     mnames.push_back(str);
    // }
    // 
    // Here's the call that tells the meta-data object that we have a mat:
    //
    // AddMaterialToMetaData(md, matname, mesh_for_mat, nmats, mnames);
    //
    //
    // Here's the way to add expressions:
    //Expression momentum_expr;
    //momentum_expr.SetName("momentum");
    //momentum_expr.SetDefinition("{u, v}");
    //momentum_expr.SetType(Expression::VectorMeshVar);
    //md->AddExpression(&momentum_expr);
    //Expression KineticEnergy_expr;
    //KineticEnergy_expr.SetName("KineticEnergy");
    //KineticEnergy_expr.SetDefinition("0.5*(momentum*momentum)/(rho*rho)");
    //KineticEnergy_expr.SetType(Expression::ScalarMeshVar);
    //md->AddExpression(&KineticEnergy_expr);
    //
}


// ****************************************************************************
//  Method: avtMOABFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: vijaysm -- generated by xml2avt
//  Creation:   Wed Jan 20 13:02:35 PST 2016
//
// ****************************************************************************

vtkDataSet *
avtMOABFileFormat::GetMesh(int domain, const char *meshname)
{

  debug1 << "avtMOABFileFormat::GetMesh from file: " << fileName << " domain: " << domain << " meshname :" << meshname  << " rank:" << rank << "\n";
  try
    {
        moab::ErrorCode merr;
        moab::Range verts, ents;

        //
        // Load the mesh file fully and populate all require information
        // related to the mesh topology and entities/sets
        //
        if (NULL==mbCore)
        {
            mbCore = new moab::Core();
            if (readOptions->FindIndex("edge") >= 0)
            {
                opt1d = readOptions->GetBool("edge");
            }
            if (readOptions->FindIndex("face") >= 0)
            {
                opt2d = readOptions->GetBool("face");
            }
            if (readOptions->FindIndex("solid") >= 0)
            {
                opt3d = readOptions->GetBool("solid");
            }
            debug1 << "avtMOABFileFormat dim options: 1d:" << (int)opt1d << " 2d:" << (int)opt2d
                    << " 3d:" << (int)opt3d <<"\n";
#ifdef PARALLEL
            string partitionMethod = readOptions->GetString("Partition:");
            string ropts="STORE_SETS_FILEIDS;PARALLEL=READ_PART;PARTITION="+partitionMethod+";";
            moab::Interface * mb = (moab::Interface *) mbCore;
            pcomm = new moab::ParallelComm(mb, VISIT_MPI_COMM);
            debug1 << "avtMOABFileFormat pcomm->comm() == VISIT_MPI_COMM " << (pcomm->comm() == VISIT_MPI_COMM) << "\n";
#else
            string ropts("STORE_SETS_FILEIDS;");
            debug1 << "avtMOABFileFormat serial MOAB\n";
#endif
            merr = mbCore->load_file(fileName, 0,ropts.c_str() );MBVIS_CHK_ERR(merr);
#ifdef PARALLEL
            // some debugging info
            moab::Range localEnts;
            pcomm->get_part_entities(localEnts);
            debug1 << " part entities : "<< localEnts.size() << "\n";
#endif
        }
        //
        // Create the unstructured mesh
        //
        vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();

        {
            // Get the list of vertices
            merr = mbCore->get_entities_by_dimension(0, 0, verts, true);MBVIS_CHK_ERR(merr);
            debug1 << "avtMOABFileFormat::GetMesh: verts.size()=" << verts.size() << "\n";
            //
            // Get the coordinates for the vertex elements.
            //
            std::vector<double> coords(verts.size()*3);
            debug1 << "avtMOABFileFormat::GetMesh: coords.size()=" << coords.size() << "\n";
            merr = mbCore->get_coords(verts, &coords[0]);MBVIS_CHK_ERR(merr);

            //
            // Populate the coordinates.  Put in 3D points with z=0 if the mesh is 2D.
            //
            vtkPoints *points  = vtkPoints::New();
            points->SetNumberOfPoints(verts.size());
            debug1 << "avtMOABFileFormat::GetMesh: points:" <<  (vtkPoints*)points << "\n";
            float *pts = (float *) points->GetVoidPointer(0);
            for (size_t i = 0; i < verts.size()*3; i++)
            {
                pts[i] = static_cast<float> (coords[i]);
            }
            coords.clear();

            ugrid->SetPoints(points);
            points->Delete();
        }

        // TODO: Do we want all entities recursively ? If we are not working on the root set
        // think about how to handle it here
        merr = mbCore->get_entities_by_handle(0, ents, false);MBVIS_CHK_ERR(merr);
        debug1 << "avtMOABFileFormat::GetMesh: ents:" <<  ents.size() << "\n";
        const moab::EntityHandle* connect;

        std::vector<vtkIdType> conn_data; // might need reordering because moab and vtk have
                                          // different order for quadratic elements
        vtkIdType vertIds[27];

        // remove the vertices and the entity sets
        moab::Range vts = ents.subset_by_type(moab::MBVERTEX);
        ents = subtract(ents, vts);
        vts = ents.subset_by_type(moab::MBENTITYSET);
        ents = subtract(ents, vts);
        debug1 << "avtMOABFileFormat::GetMesh: strip: ents:" <<  ents.size() << "\n";
        // filter by dimension
        if (opt1d) *edges = ents.subset_by_dimension(1); // edges
        if (opt2d) *faces = ents.subset_by_dimension(2); // faces
        if (opt3d) *solids = ents.subset_by_dimension(3); // solids
        *select = *edges;
        select->merge(*faces);
        select->merge(*solids);
        debug1 << "avtMOABFileFormat::GetMesh: strip: selected ents:" <<  select->size() << " psize:" << select->psize() << "\n";
        for (size_t  i = 0; i < select->size(); i++)
        {
            moab::EntityType eType = mbCore->type_from_handle((*select)[i]);

            int numnodes;

            // get all nodes now, not only corner nodes
            merr = mbCore->get_connectivity((*select)[i], connect, numnodes);MBVIS_CHK_ERR(merr);
            const moab::VtkElemType * mbToVtk = moab::VtkUtil::get_vtk_type(eType, numnodes);

            if (!mbToVtk)
            {
                debug3 << " moab type " << eType << " with " << numnodes << " is unsupported\n";
                continue; // unsupported translation from h5m to vtk; we should not have any
            }

            debug3 << " entity " << i << " mbtype " << eType << " num_nodes: " << numnodes << " vtk type " << mbToVtk->vtk_type << endl;
            if (eType!=moab::MBPOLYHEDRON)
            {
                debug3 << "  " ;
                conn_data.resize(numnodes);
                for (int j = 0; j < numnodes; ++j)
                {
                    conn_data[j] = verts.index(connect[j]);
                    debug3 << conn_data[j] << " " ;
                }
                debug3 <<endl;

                if (mbToVtk->node_order)
                {
                    debug3 << "   reorder: ";
                    for (int k = 0; k < numnodes; ++k)
                    {
                        vertIds[k] = conn_data[mbToVtk->node_order[k]];
                        debug3 << vertIds[k] << " " ;
                    }
                    debug3 <<endl;
                    ugrid->InsertNextCell(mbToVtk->vtk_type, numnodes, vertIds);
                }
                else // if there is no order change, vtk and moab agree
                {
                    ugrid->InsertNextCell(mbToVtk->vtk_type, numnodes, &conn_data[0]);
                }

            }
            else
            {
                // adapt from exodus

                vtkIdType cellarr_buf[1024]; // assume we do not have polyhedra with a lot of vertices...
                // the capacity should be OK
                int nfaces = numnodes;// for polyhedron, the connectivity contains the faces
                // the buffer will contain the
                int qq=0;
                //cellarr_buf[qq++] = nfaces;
                for (int jj = 0; jj < nfaces; jj++)
                {
                    const moab::EntityHandle* faceconn;
                    merr = mbCore->get_connectivity(connect[jj], faceconn, numnodes);MBVIS_CHK_ERR(merr);

                    cellarr_buf[qq++] = numnodes;
                    for (int kk = 0; kk < numnodes && qq<1024 ; kk++)
                    {
                        cellarr_buf[qq++] = verts.index(faceconn[kk]);
                    }
                }
                if (DebugStream::Level3())
                {
                    debug3 <<"cellar_buff: " ;
                    {
                        for (int jj=0; jj<qq; jj++)
                        {
                            debug3 << " " << cellarr_buf[jj];
                        }
                        debug3 << endl;
                    }
                }
                ugrid->InsertNextCell(VTK_POLYHEDRON, nfaces, cellarr_buf);

            }

        }

        return ugrid;
    }
    catch (moab::ErrorCode errCode)
    {
        std::string errInfo;
        mbCore->get_last_error(errInfo);
        if (!avtCallback::IssueWarning(errInfo.c_str()))
            std::cerr << errInfo << std::endl;
    }
    return 0;
}


// ****************************************************************************
//  Method: avtMOABFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: vijaysm -- generated by xml2avt
//  Creation:   Wed Jan 20 13:02:35 PST 2016
//
// ****************************************************************************

vtkDataArray *
avtMOABFileFormat::GetVar(int domain, const char *varname)
{
    if (NULL==mbCore)
      EXCEPTION1(InvalidVariableException, varname);

    try
    {

        string meshName = metadata->MeshForVar(varname);
        const avtMeshMetaData *mmd = metadata->GetMesh(meshName);
        debug1 << "avtMOABFileFormat::GetVar varname: " << varname << " meshname:" << meshName << "\n";
        // node or element dense tags?

        if (string(varname) == "ParallelPartition")
            return GetPartitionTagAsEnumScalar();
        if (string(varname) == "Materials")
            return GetMaterialTagAsEnumScalar();

        if (string(varname) == "NeumannSets")
            return GetNeumannSetsVar();

        if (string(varname) == "DirichletSets")
            return GetDirichletSetsVar();

        if (string(varname) == "GeometrySets")
            return GetGeometrySetsVar();

        bool nodeTag = (strncmp("NODE_", varname, 5 )==0);
        bool elemTag = (strncmp("ELEM_", varname, 5 )==0);

        debug1 << "avtMOABFileFormat::GetVar varname: " << varname << " is elem tag? : " << elemTag << "\n";
        string tagName;
        if (nodeTag)
        {
            tagName = string(varname+5);
        }
        else if (elemTag)
        {
            tagName = string(varname+5);
        }

        debug1 << "moab tag name:" <<tagName << "\n";
        vtkDataArray * result = 0;

        moab::Range ents;
        moab::ErrorCode merr;
        if (nodeTag)
        {
            merr = mbCore->get_entities_by_type(0, moab::MBVERTEX, ents);MBVIS_CHK_ERR(merr);
        }
        else if (elemTag)
        {
            ents = *select;
        }

        debug1 << "avtMOABFileFormat::GetVar: ents.size()=" <<  ents.size() <<  "  psize :" << ents.psize() << "\n";

        moab::Tag tag;
        merr = mbCore->tag_get_handle(tagName.c_str(),  tag);MBVIS_CHK_ERR(merr);

        // moab tag data type
        moab::DataType tag_type;
        merr = mbCore->tag_get_data_type(tag,  tag_type);MBVIS_CHK_ERR(merr);

        int tag_size;
        merr = mbCore->tag_get_length(tag, tag_size);

        int nents = (int) ents.size();
        debug1 << "avtMOABFileFormat::GetVar num entities " << nents << "type: " << tag_type << " size:" << tag_size << "\n";
        if (moab::MB_TYPE_INTEGER == tag_type)
        {
            vtkIntArray *ia = vtkIntArray::New();
            ia->SetNumberOfComponents(tag_size); // should be 1...
            ia->SetNumberOfTuples(nents);

            // use tag iterate to access direct memory
            moab::Range::iterator iter = ents.begin();
            void * data=NULL; //used for stored area
            int count =0;
            int indexInIA = 0;

            while (iter != ents.end())
            {
                merr = mbCore->tag_iterate(tag, iter, ents.end(), count, data); MBVIS_CHK_ERR(merr);
                int * ptrTag=(int*)data;
                for (int i=0; i<count; i++, ++iter, indexInIA++, ptrTag++)
                {
                    ia->SetValue(indexInIA, *ptrTag);
                }
            }
            result = ia;
        }
        if (moab::MB_TYPE_DOUBLE == tag_type)
        {
            vtkFloatArray *fa = vtkFloatArray::New();
            fa->SetNumberOfComponents(tag_size); // should be 1
            fa->SetNumberOfTuples(nents);

            // use tag iterate to access direct memory
            moab::Range::iterator iter = ents.begin();
            void * data=NULL; //used for stored area
            int count =0;
            int indexInFA = 0;

            debug1 << "avtMOABFileFormat::GetVar num entities " << nents << "type: " << tag_type << " size:" << tag_size << "\n";
            while (iter != ents.end())
            {
              debug2 << "avtMOABFileFormat::GetVar at iter: " << mbCore->id_from_handle(*iter) << " type: " << mbCore->type_from_handle(*iter)<< "\n";
              merr = mbCore->tag_iterate(tag, iter, ents.end(), count, data); MBVIS_CHK_ERR(merr);
              debug2 << "avtMOABFileFormat::GetVar count: " <<count<< "\n";
              double * ptrTag=(double*)data;
              for (int i=0; i<count; i++, iter++)
              {
                debug5 << "avtMOABFileFormat::GetVar i: " <<i << " indexInFA:" << indexInFA<< "*ptrTag " <<(*ptrTag) <<  "\n";
                for (int j=0; j<tag_size; j++, indexInFA++, ptrTag++) // really, it should be only one component
                {
                  fa->SetValue(indexInFA, (float) (*ptrTag) );

                }
              }
            }
            result = fa;
        }
        return result;
    }
    catch (moab::ErrorCode errCode)
    {
        std::string errInfo;
        mbCore->get_last_error(errInfo);
        if (!avtCallback::IssueWarning(errInfo.c_str()))
            std::cerr << errInfo << std::endl;
    }
    return 0;
}

vtkDataArray*
avtMOABFileFormat::GetPartitionTagAsEnumScalar(){
    // entities will be all of them
    // maybe we should save them all
    try {
        moab::Range ents;
        moab::ErrorCode merr = mbCore->get_entities_by_handle(0, ents, false);MBVIS_CHK_ERR(merr);
        // remove the vertices and the entity sets
        moab::Range vts = ents.subset_by_type(moab::MBVERTEX);
        ents = *select;
	if (ents.empty())
          ents = vts; 
        vtkIntArray *pparr = vtkIntArray::New();
        pparr->SetNumberOfComponents(1);
        pparr->SetNumberOfTuples(ents.size());

        int num_levels = num_parts + 1;
        const char * name = "PARALLEL_PARTITION";

        debug2 << " name:" << name << "\n";
        int * tempArr = new int [(int)ents.size()];
        for (int j=0; j<(int)ents.size(); j++)
          tempArr[j] = num_levels  ;
        // get moab tag for parallel partitions
        moab::Tag ptag;
        merr = mbCore->tag_get_handle(name, ptag);MBVIS_CHK_ERR(merr);
        // get sets of type and tag
        moab::Range tagged_sets;
        merr = mbCore->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &ptag, NULL, 1,
                      tagged_sets, moab::Interface::UNION); MBVIS_CHK_ERR(merr);

        std::vector<int>  values;
        values.resize(tagged_sets.size());
        debug2 << " num of tagged sets:" << tagged_sets.size() << "\n";
        merr = mbCore->tag_get_data(ptag, tagged_sets, &values[0]);

        // get the value for each tag, and set it to the temp array

        int k=0;
        for (moab::Range::iterator sit = tagged_sets.begin(); sit != tagged_sets.end(); sit++, k++)
        {
            moab::EntityHandle pset = *sit;
            int val= values[k]; // already read
            moab::Range setents;
            merr = mbCore->get_entities_by_handle(pset, setents);

            for (moab::Range::iterator eit=setents.begin(); eit!=setents.end(); eit++ )
            {
                moab::EntityHandle eh = *eit;
                int idx = ents.index(eh);
                if (idx>=0)
                    tempArr[idx] =  val+1;
            }
        }

        for (int j=0; j<(int)ents.size(); j++)
        {
            debug2 << " ent j " << j << " type " << mbCore->type_from_handle(ents[j]) <<
                " id " << mbCore->id_from_handle(ents[j])
                << " tempArr: " << tempArr[j] <<"\n";
            pparr->SetComponent(j, 0, tempArr[j]);
        }
        delete [] tempArr;
        return pparr;

    }
    catch (moab::ErrorCode errCode)
    {
        std::string errInfo;
        mbCore->get_last_error(errInfo);
        if (!avtCallback::IssueWarning(errInfo.c_str()))
            std::cerr << errInfo << std::endl;
    }
    return 0;
}

vtkDataArray*
avtMOABFileFormat::GetMaterialTagAsEnumScalar()
{
  // entities will be all of them
  // maybe we should save them all
  try {
      moab::Range ents;
      moab::ErrorCode merr = mbCore->get_entities_by_handle(0, ents, false);MBVIS_CHK_ERR(merr);
      // remove the vertices and the entity sets
      moab::Range vts = ents.subset_by_type(moab::MBVERTEX);
      ents = *select;
      vtkIntArray *pparr = vtkIntArray::New();
      pparr->SetNumberOfComponents(1);
      pparr->SetNumberOfTuples(ents.size());

      const char * name = "MATERIAL_SET";

      debug2 << " name:" << name << "\n";
      int * tempArr = new int [(int)ents.size()];
      for (int j=0; j<(int)ents.size(); j++)
        tempArr[j] = 0  ; // the unspecified level is 0
      // get moab tag for parallel partitions
      moab::Tag ptag;
      merr = mbCore->tag_get_handle(name, ptag);MBVIS_CHK_ERR(merr);
      // get sets of type and tag
      moab::Range tagged_sets;
      merr = mbCore->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &ptag, NULL, 1,
                    tagged_sets, moab::Interface::UNION); MBVIS_CHK_ERR(merr);

      std::vector<int>  values;
      values.resize(tagged_sets.size());
      debug2 << " num of tagged sets:" << tagged_sets.size() << "\n";
      merr = mbCore->tag_get_data(ptag, tagged_sets, &values[0]);

      // get the value for each tag, and set it to the temp array

      int k=0;
      for (moab::Range::iterator sit = tagged_sets.begin(); sit != tagged_sets.end(); sit++, k++)
      {
          moab::EntityHandle pset = *sit;
          int val= values[k]; // already read
          moab::Range setents;
          merr = mbCore->get_entities_by_handle(pset, setents, true); // material sets need true here

          for (moab::Range::iterator eit=setents.begin(); eit!=setents.end(); eit++ )
          {
              moab::EntityHandle eh = *eit;
              int idx = ents.index(eh);
              if (idx>=0)
                  tempArr[idx] =  val;
          }
      }

      for (int j=0; j<(int)ents.size(); j++)
      {
          debug2 << " ent j " << j << " type " << mbCore->type_from_handle(ents[j]) <<
              " id " << mbCore->id_from_handle(ents[j])
              << " tempArr: " << tempArr[j] <<"\n";
          pparr->SetComponent(j, 0, tempArr[j]);
      }
      delete [] tempArr;
      return pparr;

  }
  catch (moab::ErrorCode errCode)
  {
      std::string errInfo;
      mbCore->get_last_error(errInfo);
      if (!avtCallback::IssueWarning(errInfo.c_str()))
          std::cerr << errInfo << std::endl;
  }
  return 0;
}

vtkDataArray*
avtMOABFileFormat::GetNeumannSetsVar()
{
    try {
        moab::Range ents;
        moab::ErrorCode merr = mbCore->get_entities_by_handle(0, ents, false);MBVIS_CHK_ERR(merr);
        // remove the vertices and the entity sets
        moab::Range vts = ents.subset_by_type(moab::MBVERTEX);
        ents = *select;

        const int bpuc = sizeof(unsigned char)*8;
        vtkBitArray *retval = 0;
        retval = vtkBitArray::New();
        //                                      this + 1 added for universal
        retval->SetNumberOfComponents(((num_neumann+1+bpuc-1)/bpuc)*bpuc);
        retval->SetNumberOfTuples(ents.size());
        memset(retval->GetVoidPointer(0), 0, retval->GetSize()/bpuc);

        const char * name = "NEUMANN_SET";
        debug2 << " name:" << name << "\n";

        // get moab tag for neumann set
        moab::Tag ptag;
        merr = mbCore->tag_get_handle(name, ptag);MBVIS_CHK_ERR(merr);
        // get sets of type and tag
        moab::Range tagged_sets;
        merr = mbCore->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &ptag, NULL, 1,
                      tagged_sets, moab::Interface::UNION); MBVIS_CHK_ERR(merr);

        std::vector<int>  values;
        values.resize(tagged_sets.size());
        debug2 << " num of tagged sets:" << tagged_sets.size() << "\n";
        merr = mbCore->tag_get_data(ptag, tagged_sets, &values[0]);

        // get the value for each tag, and set it to the temp array

        int k=0;
        for (moab::Range::iterator sit = tagged_sets.begin(); sit != tagged_sets.end(); sit++, k++)
        {
            moab::EntityHandle pset = *sit;
            int val= values[k]; // already read
            debug2 << " neumann/dirichlet set val: " << val << "\n";
            int index = (int) std::distance(neumannsets.begin(), neumannsets.find(val));
            if (index>=neumannsets.size())
              continue; // not in values, must be an error somewhere
            moab::Range setents;
            merr = mbCore->get_entities_by_handle(pset, setents, true); // material sets, neumann sets need true here

            debug2 << " num of ents in neu set " << val <<" at index " <<
                index << " " << setents.size()<<"\n";

            for (moab::Range::iterator eit=setents.begin(); eit!=setents.end(); eit++ )
            {
                moab::EntityHandle eh = *eit;
                int idx = ents.index(eh);
                if (idx>=0)
                    retval->SetComponent(idx, index, 1);
                debug2 << " idx " << idx << " eh:" << eh << "\n";
            }
        }


        return retval;

    }
    catch (moab::ErrorCode errCode)
    {
        std::string errInfo;
        mbCore->get_last_error(errInfo);
        if (!avtCallback::IssueWarning(errInfo.c_str()))
            std::cerr << errInfo << std::endl;
    }
    return 0;
}

vtkDataArray*
avtMOABFileFormat::GetDirichletSetsVar()
{
    try {
        moab::Range ents;
        moab::ErrorCode merr = mbCore->get_entities_by_type(0, moab::MBVERTEX, ents);MBVIS_CHK_ERR(merr);

        const int bpuc = sizeof(unsigned char)*8;
        vtkBitArray *retval = 0;
        retval = vtkBitArray::New();
        // one more for universal diri
        retval->SetNumberOfComponents(((num_diri+1+bpuc-1)/bpuc)*bpuc);
        retval->SetNumberOfTuples(ents.size());
        memset(retval->GetVoidPointer(0), 0, retval->GetSize()/bpuc);

        const char * name = "DIRICHLET_SET";
        debug2 << " name:" << name << "\n";

        // get moab tag for diri set
        moab::Tag ptag;
        merr = mbCore->tag_get_handle(name, ptag);MBVIS_CHK_ERR(merr);
        // get sets of type and tag
        moab::Range tagged_sets;
        merr = mbCore->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &ptag, NULL, 1,
                      tagged_sets, moab::Interface::UNION); MBVIS_CHK_ERR(merr);

        std::vector<int>  values;
        values.resize(tagged_sets.size());
        debug2 << " num of tagged sets:" << tagged_sets.size() << "\n";
        merr = mbCore->tag_get_data(ptag, tagged_sets, &values[0]);

        int k=0;
        for (moab::Range::iterator sit = tagged_sets.begin(); sit != tagged_sets.end(); sit++, k++)
        {
            moab::EntityHandle pset = *sit;
            int val= values[k]; // already read
            debug2 << " dirichlet set val: " << val << "\n";
            int index = (int) std::distance(dirichsets.begin(), dirichsets.find(val));
            if (index>=dirichsets.size())
              continue; // not in values, must be an error somewhere
            moab::Range setents;
            merr = mbCore->get_entities_by_handle(pset, setents, true);MBVIS_CHK_ERR(merr); // material sets, neumann, diri sets need true here
            debug2 << " num of ents in diri set " << val <<" at index " <<
                index << " " << setents.size()<<"\n";

            for (moab::Range::iterator eit=setents.begin(); eit!=setents.end(); eit++ )
            {
                moab::EntityHandle eh = *eit;
                int idx = ents.index(eh);
                if (idx>=0)
                    retval->SetComponent(idx, index, 1);
                debug3 << " idx " << idx << " eh:" << eh << "\n";
            }
        }


        return retval;

    }
    catch (moab::ErrorCode errCode)
    {
        std::string errInfo;
        mbCore->get_last_error(errInfo);
        if (!avtCallback::IssueWarning(errInfo.c_str()))
            std::cerr << errInfo << std::endl;
    }
    return 0;
}

vtkDataArray*
avtMOABFileFormat::GetGeometrySetsVar()
{
   try {
        moab::Range ents;
        moab::ErrorCode merr = mbCore->get_entities_by_handle(0, ents, false);MBVIS_CHK_ERR(merr);
        // remove the vertices and the entity sets
        moab::Range vts = ents.subset_by_type(moab::MBVERTEX);
        ents = *select;

        const int bpuc = sizeof(unsigned char)*8;
        vtkBitArray *retval = 0;
        retval = vtkBitArray::New();
        // one more for universal geom
        retval->SetNumberOfComponents(((num_geom+1+bpuc-1)/bpuc)*bpuc);
        retval->SetNumberOfTuples(ents.size());
        memset(retval->GetVoidPointer(0), 0, retval->GetSize()/bpuc);

        const char * name = "GEOM_DIMENSION";
        debug2 << " name:" << name << "\n";

        // get moab tag for geom set
        moab::Tag ptag;
        merr = mbCore->tag_get_handle(name, ptag);MBVIS_CHK_ERR(merr);
        // get sets of type and tag
        moab::Range tagged_sets;
        merr = mbCore->get_entities_by_type_and_tag(0, moab::MBENTITYSET, &ptag, NULL, 1,
                      tagged_sets, moab::Interface::UNION); MBVIS_CHK_ERR(merr);

        // now get the file ids on all sets
        moab::Tag setFileIdTag;
        merr = mbCore->tag_get_handle("__FILE_ID_FOR_SETS", setFileIdTag); MBVIS_CHK_ERR(merr);

        std::vector<long> setsFileIds(tagged_sets.size());
        merr = mbCore->tag_get_data(setFileIdTag, tagged_sets, &setsFileIds[0]); MBVIS_CHK_ERR(merr);

        // subtract the set start file id
        debug2 << "sets ids: \n";
        for (size_t i=0; i<setsFileIds.size(); i++)
        {
            setsFileIds[i] = setsFileIds[i] - file_descriptor->sets.start_id +1 ;
            debug2 << "i, setsFileIds[i] = " << i << " " << setsFileIds[i] << "\n";
        }

        std::vector<int> values;
        values.resize(tagged_sets.size());
        debug2 << " num of tagged sets:" << tagged_sets.size() << "\n";
        merr = mbCore->tag_get_data(ptag, tagged_sets, &values[0]); MBVIS_CHK_ERR(merr);

        int k=0;
        for (moab::Range::iterator sit = tagged_sets.begin(); sit != tagged_sets.end(); sit++, k++)
        {
            moab::EntityHandle pset = *sit;
            int fileId = (int) setsFileIds[k];
            int index = std::lower_bound( &(file_descriptor->defTagsEntSets[4][0]),
                &(file_descriptor->defTagsEntSets[4][0]) + num_geom, fileId) - &(file_descriptor->defTagsEntSets[4][0]);
            if (index >= num_geom)
              continue;
            int val= values[k]; // already read
            debug2 << " geometry set val: " << val << " index " << index<< " fileid " << setsFileIds[k]   << "\n";

            moab::Range setents;
            merr = mbCore->get_entities_by_handle(pset, setents, true);MBVIS_CHK_ERR(merr); // material sets, neumann, diri sets need true here
            debug2 << " num of ents in geom set " << val <<" at index " <<
                index << " " << setents.size()<<"\n";

            for (moab::Range::iterator eit=setents.begin(); eit!=setents.end(); eit++ )
            {
                moab::EntityHandle eh = *eit;
                int idx = ents.index(eh);
                if (idx>=0)
                    retval->SetComponent(idx, index, 1);
                debug2 << " idx " << idx << " eh:" << eh << "\n";
            }
        }


        return retval;

    }
    catch (moab::ErrorCode errCode)
    {
        std::string errInfo;
        mbCore->get_last_error(errInfo);
        if (!avtCallback::IssueWarning(errInfo.c_str()))
            std::cerr << errInfo << std::endl;
    }
    return 0;
}
// ****************************************************************************
//  Method: avtMOABFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: vijaysm -- generated by xml2avt
//  Creation:   Wed Jan 20 13:02:35 PST 2016
//
// ****************************************************************************
vtkDataArray *
avtMOABFileFormat::GetVectorVar(int domain, const char *varname)
{
  if (NULL==mbCore)
    EXCEPTION1(InvalidVariableException, varname);

  try {

    string meshName = metadata->MeshForVar(varname);
    const avtMeshMetaData *mmd = metadata->GetMesh(meshName);
    debug1 << "avtMOABFileFormat::GetVectorVar varname: " << varname << " meshname:"
        << meshName << "\n";
    // node or element dense tags?

    bool nodeTag = (strncmp("NODE_", varname, 5) == 0);
    bool elemTag = (strncmp("ELEM_", varname, 5) == 0);

    debug1 << "avtMOABFileFormat::GetVectorVar varname: " << varname
        << " is elem tag? : " << elemTag << "\n";
    string tagName;
    if (nodeTag) {
      tagName = string(varname + 5);
    } else if (elemTag) {
      tagName = string(varname + 5);
    }

    debug1 << "moab tag name:" << tagName << "\n";
    vtkDataArray * result = 0;

    moab::Range ents;
    moab::ErrorCode merr;
    if (nodeTag) {
      merr = mbCore->get_entities_by_type(0, moab::MBVERTEX, ents);
      MBVIS_CHK_ERR(merr);
    } else if (elemTag) {
      ents = *select;
    }

    moab::Tag tag;
    merr = mbCore->tag_get_handle(tagName.c_str(), tag);
    MBVIS_CHK_ERR(merr);

    // moab tag data type
    moab::DataType tag_type;
    merr = mbCore->tag_get_data_type(tag, tag_type);
    MBVIS_CHK_ERR(merr);

    int tag_size;
    merr = mbCore->tag_get_length(tag, tag_size);

    int nents = (int) ents.size();
    debug1 << "avtMOABFileFormat::GetVectorVar num entities " << nents << "type: "
        << tag_type << " size:" << tag_size << "\n";
    if (moab::MB_TYPE_INTEGER == tag_type) {
      vtkIntArray *ia = vtkIntArray::New();
      ia->SetNumberOfComponents(tag_size); // should be > 1
      ia->SetNumberOfTuples(nents);

      // use tag iterate to access direct memory
      moab::Range::iterator iter = ents.begin();
      void * data = NULL; //used for stored area
      int count = 0;
      int indexInIA = 0;

      while (iter != ents.end()) {
        merr = mbCore->tag_iterate(tag, iter, ents.end(), count, data);
        MBVIS_CHK_ERR(merr);
        int * ptrTag = (int*) data;
        for (int i = 0; i < count; i++, iter++)
          for (int j = 0; j < tag_size; j++, indexInIA++, ptrTag++) {
            ia->SetValue(indexInIA, (float) (*ptrTag));
          }
      }
      result = ia;
    }
    if (moab::MB_TYPE_DOUBLE == tag_type) {
      vtkFloatArray *fa = vtkFloatArray::New();
      fa->SetNumberOfComponents(tag_size); // could be more than 1 !!
      fa->SetNumberOfTuples(nents);

      // use tag iterate to access direct memory
      moab::Range::iterator iter = ents.begin();
      void * data = NULL; //used for stored area
      int count = 0;
      int indexInFA = 0;

      debug1 << "avtMOABFileFormat::GetVectorVar num entities " << nents << "type: "
          << tag_type << " size:" << tag_size << "\n";
      while (iter != ents.end()) {
        debug2 << "avtMOABFileFormat::GetVectorVar at iter: "
            << mbCore->id_from_handle(*iter) << " type: "
            << mbCore->type_from_handle(*iter) << "\n";
        merr = mbCore->tag_iterate(tag, iter, ents.end(), count, data);
        MBVIS_CHK_ERR(merr);
        debug2 << "avtMOABFileFormat::GetVectorVar count: " << count << "\n";
        double * ptrTag = (double*) data;
        for (int i = 0; i < count; i++, iter++) {
          debug5 << "avtMOABFileFormat::GetVectorVar i: " << i << " indexInFA:"
              << indexInFA << "*ptrTag " << (*ptrTag) << "\n";
          for (int j = 0; j < tag_size; j++, indexInFA++, ptrTag++) {
            fa->SetValue(indexInFA, (float) (*ptrTag));

          }
        }
      }
      result = fa;
    }
    return result;
  } catch (moab::ErrorCode errCode) {
    std::string errInfo;
    mbCore->get_last_error(errInfo);
    if (!avtCallback::IssueWarning(errInfo.c_str()))
      std::cerr << errInfo << std::endl;
  }
  return 0;
}
